The  above  pictures  show  the  analytical  and  computational  aspects  of  turbulence  and  combustion.  These  figures  show  an  un¬ 
stable  detonation  wave  in  two  distinct  states.  The  top  figure  snows  the  temperature  profile  of  a  Z-N-D  detonation  wave  with  qo 
(the  heat  release  parameter)  of  50,  E*  (^e  activation  energy)  of  10,  and  f  (the  overdriver  of  the  detonation)  of  1.2.  Note  the 
generatirxt  of  regularly  paired  eddies.  The  bottom  picture  is  produced  under  the  same  parameter  values  except  for  E*  that 
has  a  value  of  50.  This  picture  shows  the  generation  of  strong  irregular  cells  and  fully  developed  2-D  turbulence. 

This  pioneering  research  of  Professor  Andrew  Majda  and  co-workers  at  Princeton  University  ivas  sponsored  by  the  Office  of 
Nava  Researm.  On  page  20,  Professor  Majda  is  featured  in  "Profiles  in  Science. " 

ONR  /ms  been  working  in  this  field  of  research  for  30  years.  The  discovery  in  the  1960^  under  ONR  support  of  coherence  in 
turbulence  revotutiorvzed  the  view  of  turbulent  flow.  Before  that  time  the  statistical  description  of  turbulence  emphasized  its 
randomness  and  brought  about  time-averaged  turbulence  models.  ONR-sponsored  research  identified  coherent  structures  in 
both  free  and  bounded  shear  flows,  and  subsequent  prograns  investigated  the  role  of  these  structures  in  the  production  ana 
dissipation  oT  energy  in  the  flow. 
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Research  on  Multiple 
Valued  Logic  At  The 
Naval  Postgraduate 
School 

Jon  T.  Butler 

Department  of  Electrical  and  Computer  Engineering 
Naval  Postgraduate  School 


Abstract 

The  Navy’s  need  for  faster,  smaller  computers  has  in¬ 
spired  research  on  the  use  of  more  than  two  levels  of  logic  in 
computer  circuits.  This  paper  focusses  on  work  at  NPS  in 
multiple-valued  logic.  It  discusses  the  development  of  multi¬ 
ple-valued  programmable  logic  arrays  (MVL-PLA),  as  well 
as  computer-aided  design  tools  for  MVL-PLA’s. 

Introduction 

The  effectiveness  of  Navy  systems  is  dependent  on  com¬ 
puters.  Computers  are  used  to  predict  weather,  to  track  re¬ 
sources,  to  guide  missiles,  to  simulate  combat,  and  to  process 
the  words  and  numbers  necessary  for  the  general  conduct  of 
Navy  business.  Their  usefulness  depends  on  computing 
powCT,  speed,  and  sue.  As  with  other  critical  resources,  the 
Navy  cannot  depend  on  the  private  sector  as  the  exclusive 
provider  of  the  technology:  it  must  be  an  active  participant. 
This  article  focuses  on  computer  circuits,  and  specifically  on 
multiple-valued  circuits,  whose  use  improves  computing 
power,  speed,  and  size.  We  first  consider  binary  circuits.  Then, 


we  discuss  multiple-valued  circuits.  This  serves  as  introduc¬ 
tion  to  the  main  discussion,  research  at  the  Naval  Postgraduate 
School  on  this  important  subject. 

Disadvantages  of  Binary 

Current  computers  are  binary;  they  use  two  values  of 
logic.  For  convenience,  we  represent  these  as  0  and  1 .  But  0 
and  1  are  simply  symbols  for  some  quantity,  such  as  voltage, 
current,  or  charge.  For  example,  in  the  most  widely  used 
circuits,  0  represents  0  volts  and  1  represents  5  volts. 

Binary  is  commonly  used  for  historical  reasons.  Comput¬ 
ers  in  the  1940’s  used  relays,  which  had  two  stable  stales,  open 
and  closed.  Tubes  nrd  transistors  have  two  stable  states,  satu¬ 
ration  (conducting)  and  cutoff  (nonconducting).  There  are  dis¬ 
advantages  to  binary,  however.  One  is  evident  by  the  number  of 
bits  needed  to  represent  a  number.  For  example,  the  decimal 
number  33  corresponds  to  10001  in  binary.  That  is,  while  five 
symbols  are  needed  to  represent  33  in  binary,  only  two  are  needed 
in  decimal.  Unfortunately,  the  difference  increases  with  the  size 
of  the  number.  This  inherent  complexity  is  hidden  from  most 
users  by  hardware  (e.g.,  binary  to  decimal  converters  in  digital 
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watches)  and  by  software  (e.g.,  high  level  languages).  How¬ 
ever,  certain  disadvantages  cannot  be  so  simply  hidden.  Real 
penalties  exist  in  both  compactness  and  speed. 

Compactness 

The  complexity  of  the  circuits  that  perform  addition, 
multiplication,  and  other  arithmetic  operations  depends  on  the 
number  of  symbols  used.  For  commonly  used  number  sys¬ 
tems,  this  complexity  increases  with  the  number  of  symbols. 
As  a  result,  so  does  the  chip  area  needed  for  devices.  However, 
the  most  significant  praalty  is  not  with  the  devices,  but  the 
interconnect  between  devices.  In  VLSI  circuits,  about  70%  of 
the  chip  area  is  used  for  interconnect,  20%  for  insulation,  and 
10%  for  devices.  Interconnect  corresponds  to  the  chip  area 
used  to  carry  information  around  the  circuit.  Insulation  is 
needed  to  separate  devices  and  interconnect.  With  so  much 
area  used  to  connect  devices,  less  area  is  available  for  logic. 

The  inefTicient  use  of  interconnect  by  binary  also  occurs 
outsitte  the  chip.  Signals  going  into  and  out  of  a  VLSI  chip 
must  use  pins  connecting  the  chip  to  the  circuit  board.  This  is 
called  the  pinout  problem.  That  is,  while  significant  reductions 
have  occurred  in  the  size  of  logic  devices,  there  has  been  no 
comparable  reduction  in  the  size  of  integrated  circuit  pin 
connections.  Strength  and  reliability  dictate  a  (relatively  laige) 
minimum  pin  connection  size.  As  a  result,  there  is  a  need  to 
better  use  existing  pins. 

Speed 

A  limit  on  the  speed  of  arithmetic  circuits  is  the  carry  (or 
borrow)  betwera  digits.  For  example,  in  binary  addition,  the 
most  significant  sum  bit  depends  not  only  on  the  most  signif¬ 
icant  bits  of  the  two  numbers  being  added,  but  also  on  all  lower 
bits  (contrast  this  to  the  least  significant  sum  bit  which  does 
not  depend  on  higher  order  bits).  The  dependence  occurs 
through  the  carry  between  bits.  Thus,  computation  of  the  most 
significant  sum  bit  must  wait  until  the  carries  out  of  all  lower 
order  bits  are  computed.  Addition,  when  done  in  this  way, 
cannot  be  performed  at  high  speed. 

However,  there  are  multiple-valued  number  systems 
where  this  dependence  does  not  exist  An  example  is  the 
residue  number  system.  In  this  system,  operations  that  form 
the  sum,  difference,  or  product  occur  only  at  each  digit  inde¬ 
pendent  of  all  other  digits.  Because  of  this  independence, 
arithmetic  operations  are  fast  in  the  residue  number  system. 

Multiple-Valued  Logic 

All  of  the  disadvantages  discussed  above  can  be  alleviated 
by  using  more  than  two  levels  of  logic.  The  problem  of  a  large 
number  of  bits  needed  to  rqrtesent  a  number  decreases  as  the 
number  of  logic  values  increases.  Interconnect,  both  internal  and 
external  to  the  VLSI  chip,  is  mote  efficiently  used  in  multiple- 


Figure  1. 

Photomicrograph  of  a  Basic  Charge  Coupled  Device  Program¬ 
mable  Logic  Array  (Courtesy  of  the  University  of  Twente.©IEEE 
1988). 


valued  logic.  The  carry  prc^gation  problem  disq^pears  with 
the  choice  of  multiple-valued  number  system. 

However,  there  remains  the  question  of  implementing  a 
multiple- valued  system.  Modem  computing  devices  naturally 
allow  more  than  two  levels  of  logic.  For  example,  in  charge- 
coupled  devices  (CCD),  logic  levels  ate  represented  as  various 
quantities  of  charge^.  Resonant-tunnelling  devices^  allow  four 
or  more  levels  of  voltage.  Fuzzy  logic  devices  offer  large 
numbers  of  logic  levels;  so  many,  that  they  are  modelled  as 
inTtnite;  present  circuit  implementations  exceed  IS  levels'^. 
Bio-molecular  devices,  where  information  processing  occurs 
at  the  molecule  level^,  operate  on  an  enormous  number  of  logic 
levels,  conservatively  estimated  to  be  over  10,000. 

Multiple-Valued  Programmable  Logic  Arrays 

The  design  of  logical  circuits  is  a  complex  process.  The 
vast  amount  of  circuitry  needed  for  even  a  simple  operation 
like  addition  can  be  daunting.  Because  of  this,  a  way  is  needed 
to  organize  circuits.  The  programmable  logic  array  or  PLA  is 
a  circuit  sUDcture  designed  to  handle  this  complexity.  PLA's 
are  uniform  and  thus  easily  replicated  in  VLSI.  They  can  be 
programmed,  and  therefore  adapted  to  realize  a  given  design 
specification.  Because  of  this  characteristic,  our  research  at  the 
Naval  Postgraduate  School  has  concentrated  on  PLA  design. 

In  a  cooperative  venture  with  the  University  of  Twente  1.. 
Enschede,  The  Netherlands,  we  have  developed  a  PLA  circuit 
for  charge-coupled  device  (CCD)  technology  (See  Figure  1 
from^.  In  CCD,  values  are  stored  in  cells  (capacitors)  as 
charge.  A  logic  0,  1,  2,  and  3  is  stored  as  0,  1,000,000, 
2,000,000,  and  3,000,000  electrons.  Computation  occurs  as 
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charge  moves  among  cells  combining  and  breaking  up.  An 
advantage  of  CCD  is  cornpacmess,  more  so  than  any  other 
VLSI  technology.  While  slower  than  the  common  CMOS 
technology,  it  is  faster  than  disk  and  stands  as  a  replacement 
for  disk.  The  use  of  multiple-valued  logic  in  CCD  increases 
its  density  significantly.  CCD  is  used  as  a  light  receptor,  and 
the  advent  of  logic  circuits  in  CCD  means  that  image  process¬ 
ing  can  be  done  on  the  same  chip  in  which  the  image  is 
received.  With  large  numbers  of  logic  levels,  there  is  no  need 
for  encoders  between  light  receptors  and  the  logic  used  to 
process  and  store  the  image.  CCD  circuits  with  many  logic 
levels  have  been  fabricated.  For  example,  Hitachi  has  suc¬ 
ceeded  in  building  a  16  level  CCD  memory. 

While  the  uniform  structure  of  a  PL  A  helps  one  to  handle 
the  high  complexity  of  logic  circuits,  there  is  the  problem  of 
designing  large  PLA’s.  To  make  multiple- valued  PLA  design 
feasible  for  large  systems,  we  initiated  a  program  to  develop 
computer-aided  design  tools. 

Computer-Aided  Design  for 
Multiple-Valued  Logic 

We  began  this  program  with  a  study  of  heuristic  tech¬ 
niques  for  finding  minimal  PLA’s,  PLA’s  with  the  smallest 
structure.  A  heuristic  will  produce  a  design,  although  not 
necessarily  a  minimal  one.  In  1987  at  the  beginning  of  this 
study,  there  were  three  heuristic  techniques.  At  that  time,  there 
was  no  systematic  analysis  of  their  relative  merits;  each  ex¬ 
isted  separately.  In  an  effort  to  gain  an  understanding  of  the 
multiple-valued  PLA  design  process,  we  compared  the  three 
heuristics  on  the  same  functions  (specifications).  One  interest¬ 
ing  result‘d  was  that  heuristics  tended  to  perform  better  on 
different  classes  of  functions,  and  there  was  a  surprising 
advantage  to  applying  all  three  and  choosing  the  most  compact 
realization. 

We  also  provided  a  partial  answer  to  an  important  ques¬ 
tion;  How  well  do  heuristics  perform  compared  to  the  best  that 
one  could  possibly  do?  To  answer  this,  we  developed  the  Hrst 
practical  approach  to  finding  provably  minimal  solutions;  an 
efHcient  search  technique  that  is  much  faster  than  exhaustive 
search,  but  still  guaranteeing  a  minimal  solution.  Indeed, 
without  the  improved  search  technique,  we  would  have  been 
unable  to  achieve  the  minimal  solution  because  of  the  enor¬ 
mous  computation  time  required  by  exhaustive  search.  This 
comparison  provided  much  needed  insight  into  heuristic  meth¬ 
ods  and  provided  realism  in  our  expectation  of  the  quality  of 
solutions  provided  by  heuristic  methods. 

It  is  imptntant  to  state  that  these  experiments  deal  with 
small  functions;  for  even  medium  size  functions  this  experi¬ 
ment  would  have  been  impossible.  For  larger  functions,  the 
only  ^Jproach  is  heuristic.  Indeed,  our  research  shows  that 
often  a  minimal  realization  is  not  produced.  Although  a  de¬ 
signer  will  not  know  how  close  his/her  design  is  to  the  opti- 


Figure  2. 


LTJG  Ugur  Ozkan  Using  HAMLET. 


mum,  our  study  indicates  that  deviations  are  probably  less  than 
10%  away.  Ours  was  the  first  systematic  study  of  the  problem. 

Our  experience  from  this  analysis  showed  the  need  for  a 
tool  to  analyze  PLA  design  heuristics.  In  1988,  we  launched 
an  extensive  design  project  for  such  a  tool.  The  result  was 
HAMLET  (Heuristic  Analyzer  for  Multiple- Valued  Logic  Ex¬ 
pression  Translation),  the  first  practical  computer-aided 
(CAD)  tool  for  multiple-valued  logic.  HAMLET'*  has  the 
ability  to  analyze  heuristics  specified  by  the  user,  using  either 
individual  functions  or  randomly  generated  functions.  This  is 
a  particularly  useful  feature.  Instead  of  writing  a  separate 
program  to  generate  test  inputs,  this  is  done  automatically  in 
HAMLET.  Of  special  significance  is  HAMLET’s  ability  to 
collect  and  automatically  display  heuristics.  Tasks  that  nor¬ 
mally  took  man-days  can  now  be  literally  done  in  keystrokes. 
HAMLET  can  produce  ensembles  of  test  cases  whose  statis¬ 
tical  properties  can  be  controlled.  The  search  algorithm  used 
in  reference  13  was  also  incorporated  into  HAMLET.  We  are 
now  able  to  use  the  search  to  find  solutions  that  were  impossi¬ 
ble  to  obtain  with  any  previous  method.  The  search  can  be 
controlled.  At  one  extreme,  it  can  go  directly  to  a  PLA  realiza¬ 
tion  of  a  given  function.  At  the  other,  it  can  perform  a  complete 
search,  producing  a  known  minimal  solution.  In  effect,  this  is 
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like  a  dial,  one  end  of  which  corresponds  to  fast  but  poor 
solutions,  while  the  other  end  corresponds  to  slow  but  minimal 
solutions.  HAMLET  can  also  do  design,  producing  the  layout 
of  a  PLA,  given  a  specification.  Figure  3  shows  one  PLA 
layout  produced  by  HAMLET.  HAMLET  is  designed  to  have 
a  long  lifetime.  It  is  easily  modified;  indeed  an  improved 
heuristic  has  been  added*'’,  as  well  as  an  entirely  new  ap¬ 
proach^  (see  below).  Also,  HAMLET  has  been  used  by  re¬ 
searchers  in  at  least  four  countries. 

Simulated  Annealing 

When  metal  (h:  glass  is  heated  and  then  slowly  cooled,  it 
is  said  to  undergo  annealing.  The  final  cooled  state  is  similar 
to  a  crystal  and  represents  a  low  energy  state.  If  it  is  cooled  at 
a  high  rate,  atoms  have  less  of  a  chance  to  align  and  a  higher 
energy  state  occurs. 

In  19S3,  Metropolis*'  proposed  a  method  of  optimizing 
complex  problems  that  mimics  this  process.  It  is  called  simu¬ 
lated  annealing.  We  have  ftmiulated  simulated  annealing  for 
the  PLA  minimization  problem  as  follows.  One  can  satisfy  a 
given  PLA  specification  by  producing  a  set  of  implicants.  Each 
implicant  represents  a  part  of  the  whole  PLA.  Solving  the  PLA 
minimization  problem,  therefore,  is  equivalent  to  Ending  the 
fewest  implicants  that  satisfies  the  speciEcation.  Indeed,  there 
is  a  direct  relationship  between  the  number  of  implicants  and 
the  chip  area  occupied.  Thus,  reducing  the  number  of  im¬ 
plicants  by  half,  for  example,  results  in  almost  a  one-half 
reduction  in  chip  area.  Implicants  can  be  manipulated.  For 
example,  one  implicant  can  be  divided  into  two,  pairs  some¬ 
times  can  be  combined  into  a  single  implicant,  and  pairs  can 
sometimes  be  changed  into  another  pair,  or  a  triple,  etc.. 
Ideally,  one  would  like  to  combine  a  pair  of  implicants  into 
one,  since  this  reduces  PLA  chip  area.  However,  in  a  given  set 
of  implicants,  it  is  quite  possible  that  no  pairs  combine.  Such 
a  set  is  said  to  be  a  local  minimum.  However,  by  replacing  a 
pair  by  another  pair,  or  by  a  triple,  etc.,  it  is  possible  to  create 
a  new  set  of  implicants,  where  now  many  combine,  leading  to 
a  smaller  set  'The  descriptor  “local”  means  there  are  no  im¬ 
plicant  sets  nearby  that  contain  fewer  elements.  It  is  not  a 
“global”  minima  if  there  is  another  set  of  implicants  (some- 
wh^)  that  has  fewer  elements. 

The  derivation  of  a  minimal  set  is  a  difficult  problem.  A 
heuristic  generates  a  solution  to  a  problem ,  but  not  necessarily 
a  minimal  one.  The  PLA  minimization  problem  belongs  to  the 
set  of  NP-complete  problems,  which  means  that  the  best 
known  algorithm  for  Ending  the  minimal  solution  increases 
exponentially  in  complexity  as  the  problem  size  increases.  It 
is  for  this  reason  that  there  is  so  much  interest  in  heuristics. 
Simulated  annealing  has  the  beneEt  that,  under  certain  condi¬ 
tions,  a  minimal  solution  will  be  found  with  a  probability  that 
approaches  100%  as  the  process  continues.  This  is  unlike  any 
other  known  heuristic  for  PLA  minimization,  where,  for  cer- 


Figure  3. 

PLA  Layout  Produced  by  HAMLET. 


tain  PLA  specifications,  the  probability  of  finding  a  minimal 
solution  is  0%. 

In  eEect,  simulated  annealing  is  a  search  through  a  space 
of  solutions.  'The  basic  computation  is  a  move.  That  is,  given 
one  solution  to  the  PLA  specification,  another  solution  is 
obtained  by  making  a  move.  Each  move  begins  by  first  select¬ 
ing  a  pair  of  implicants.  If  the  two  implicants  can  combined, 
they  are,  completing  the  move.  If  the  implicants  cannot  be 
combined,  a  replacement  is  considered  with  the  fewest  number 
of  equivalent  implicants.  In  efiect,  a  weighted  coin  is  tossed. 
If  the  coin  lands  with  heads  up,  the  replacement  is  made, 
completing  the  move.  Otherwise,  it  is  not  made,  and  another 
pair  is  selected.  The  analogy  with  annealing  is  this.  At  each 
temperature,  the  weighting  of  the  coin  is  the  same  for  some 
fixed  number  of  moves.  In  our  programs,  we  allow  20,000 
moves  at  each  temperature.  At  high  temperatures,  the  weight 
is  such  that  heads  almost  always  comes  up,  which  means  that 
all  moves  are  accepted.  However,  at  low  temperatures,  the 
weight  is  shifted  to  tails,  so  that  fewer  moves  are  accepted 
which  increase  the  number  of  implicants.  At  any  temperature, 
a  move  that  reduces  the  number  of  implicants  is  accepted. 
Thus,  at  low  temperatures,  the  mixture  of  generated  moves  is 
predominantly  PLA  complexity-decreasing  moves.  The  exis¬ 
tence  of  complexity-increasing  moves  at  all  temperatures  is 
imptMtant;  it  allows  the  system  to  escape  from  local  minima. 
This  is  analogous  to  annealing,  in  which  atoms  at  high  tem¬ 
perature  move  about  freely.  At  lower  temperatures  the  freedom 
is  restricted,  as  they  settle  into  a  low  energy  state. 

Figure  4  shows  how  simulated  annealing  performs  on  a 
multiple-valued  function  that  begins  with  96  implicants.  This 
solution  was  obtained  by  applying  a  heuristic  to  some  given 
PLA  specification.  The  axis  ladled  “number  of  product 
terms”  measures  the  quality  of  the  solution  to  the  PLA  speci¬ 
fication;  the  fewer  the  better.  The  axis  marked  temperature 
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Simulated  Annealing  Applied  to  Multiple-Valued  PLA  Minimization. 
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MULTIPLE-VALUED  PLA  MINIMIZATION 
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20.  TO  40. 
0.  TO  20. 
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RESHAPE 


measures  the  probability  with  which  implicant  increasing 
moves  are  accepted.  At  higher  temperatures  (toward  the  front), 
almost  all  implicant-increasing  moves  are  accepted.  The  ver¬ 
tical  axis  represents  the  number  of  times  a  solution  is  found  at 
the  temperature  and  the  number  of  implicants  on  the  x-y  plane. 
Here,  color  encodes  the  vertical  deviation.  In  effect,  the  verti¬ 
cal  axis  shows  how  the  simulated  annealing  program  is  doing 
with  respect  to  the  number  of  product  terms.  Specifically,  at 
the  highest  temperature  (the  slice  closest  to  the  front),  one  can 
see  how  the  simulated  annealing  progresses  from  the  initial 
solution  of  96  implicants  up  to  as  many  as  2 16  implicants.  This 
is  described  as  melting.  As  the  temperature  decreases  from  this 
point,  there  is  steady  progress  towards  PLA  specifications  with 
fewer  implicants.  In  this  application  of  simulated  annealing,  a 
solution  of  88  implicants  is  achieved.  The  behavior  of  the 
implicants  is  similar  to  the  behavior  of  individual  atoms  as 
metal  or  glass  anneals;  as  the  temperature  decreases,  the  form 


of  the  final  material  takes  shape  as  individual  atoms  occupy  a 
site  in  the  structure  corresponding  to  low  energy. 

A  comparison  of  simulated  annealing  with  other  PLA 
minimization  heuristics  shows  that  it  is  superior  on  an  average 
case  basis.  That  is,  using  HAMLET,  random  functions  were 
generated,  heuristics  were  applied  to  all,  and  the  results  com¬ 
pared.  Simulated  annealing  surpasses  the  all  known  heuristics 
on  the  solutions  obtained.  Although,  more  time  is  needed, 
simulated  annealing  has  the  benefit  that  its  performance  can 
be  controlled.  That  is,  with  larger  computation  times,  one  can 
achieve  a  more  compact  realization. 

Conclusions 

The  Navy  has  a  tradition  of  commiunent  to  computing,  as 
epitomized  by  the  contributions  of  Rear  Adm.  Grace  Murray 
Hopper.  Indeed,  in  the  past  SO  years,  significant  improvement 
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has  occurred  in  Navy  systems  because  of  the  computer.  It  is 
clear  that  the  next  SO  years  will  bring  even  further  improve¬ 
ment  The  promise  of  compact,  high  speed  computers  has 
inspired  our  work  on  multiple-valued  logic. 

Multiple- valued  logic  offers  a  means  to  overcome  signif¬ 
icant  disadvantages  of  binary.  It  makes  more  efficient  use  of 
chip  area  by  encoding  more  information  in  the  wires  that 
connect  devices,  and  it  allows  the  use  of  high-speed  carry-less 
operations.  Our  work  at  the  Naval  Postgraduate  School  has 
b^n  focused  on  develqiment  of  techniques  (PLA’s)  and  the 
tools  (CAD)  needed  to  make  multiple-vsdued  logic  a  reality. 

Our  current  effort  focuses  on  expanding  simulated  an¬ 
nealing.  This  includes  a  promising  tqiproach  to  high-speed 
computations;  indeed,  initial  indications  show  a  modification 
of  the  conventional  simulated  annealing  paradigm  can  pro¬ 
duce  computation  speeds  higher  than  any  known  heuristic, 
while  maintaining  almost  the  same  capability  to  minimize  chip 
area.  We  are  also  pursuing  the  use  of  parallel  computers  to 
design  PLA’s'^,  including  simulated  annealing. 
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Subsonic  and  Supersonic 
Mixing  and  Combustion 
Enhancement 

Gabriel  D.  Roy 
Office  of  Naval  Research 


Abstract 

A  Hve  year  research  program  was  sponsored  by  the  US 
E>epartment  of  the  Navy  to  obtain  the  science  base  for  the 
conceptual  evaluation  of  higher  speed,  long  range  supersonic 
«.cmbustion  ramjets.  This  program  was  focused  to  determine 
the  turbulent  structure,  size  scales,  intensities  and  rates  of  heat 
and  mass  transfer  of  compressible  subsonic  and  supersonic 
shear  layer  mixing  flows,  with  and  without  combustion.  Al¬ 
though  the  program  addressed  the  issues  relating  to  the  com¬ 
bustion  of  solid  fuels  in  supersonic  combustion  as  well,  the 
present  pap^  will  focus  on  supersonic  shear  layer  mixing  and 
combustion  enhancement  The  studies  showed  that  heat  re¬ 
lease  tends  to  inhibit  the  mixing  of  shear  layers.  However, 
three-dimensional  excitation  at  subharmonics  of  the  shear 
layer,  and  streamwise  vorticity  stirring  by  lobbed  splitter 
plates  have  been  found  to  increase  mixing  and  combustion. 
The  supersonic  shear  layers  are  found  to  be  stable,  and  their 
growth  deceases,  as  the  convective  Mach  number  (Me)  in¬ 
creases.  Vortex  generators  increase  the  growth  rate  of  super¬ 
sonic  shear  layers  by  about  30%,  whereas  shock  impingement 
has  very  little  effect.  Non-axisymmetric  nozzles  and  inlets 
have  shown  to  enhance  mixing  and  combustion  of  supersonic 
coaxial  flows. 


Introduction 

As  longer  ranges  are  desired  for  air  launched  tactical 
missiles,  renewed  interest  has  emerged  on  air-breathing  pro¬ 
pulsion  systems.  Aircraft  and  ship-imposed  constraints  on 
length,  diameter  and  weight  of  weanons  require  compact 
ramjets  and  scramjets.  However,  scrainjet  combustors  have 
been  suffering  decreased  mixing  associated  with  supersonic 
compressible  shear  layers,  and  the  associated  difficulty  in 
achieving  complete  combustion  in  the  limited  residence  times 
available.  To  wldress  the  fundamental  supersonic  mixing  and 
combustion  processes,  the  Office  of  Naval  Research  (ONR) 
initiated  a  focused  five  year  research  program.  The  experimen¬ 
tal,  analytical  and  numerical  research  program  addressed  the 
conceptual  evaluation  of  higher  speed,  longer  range  super¬ 
sonic  combustion  ramjets  and  the  potential  for  several  fold 
enhancement  of  the  radiant  output  of  high  performance  air¬ 
craft-launched  pyrotechnic  infrared  countermeasures.  In  par¬ 
ticular,  the  initiative  investigated  the  temporal  and  spatial 
scales,  and  rates  of  subsonic  and  supersonic  mixing  and  com¬ 
bustion  processes,  the  supersonic  flame  structures  in  mixing 
layers  adjacent  to  bounding  solid  surfaces,  and  the  character¬ 
istics  of  the  gasification  and  particle  ejection  processes  from 
high  energy  polymeric  solid  fuels  and  their  boron/metal  filled 
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Figur*  1. 

Product  Concentration  without  and  with  Energy  Release 


as  in  those  compressible  subsonic  shear  layer  mixing 
flows? 

3.  Can  the  supersonic  mixing  layers  be  manipulated  by  low 
level  forcing  at  subharmonics  of  the  natural  frequency 
of  the  mixing  layers?  Are  there  other  methods  of  con¬ 
trol? 

4.  Can  the  compressible  subsonic  and  supersonic  mixing 
layers  and  flow  structures  be  influenced  by  combustion 
generated  heat  release? 

This  paper  is  limited  to  the  research  conducted  under  this 
special  focus  program,  and  complimentary  research  under  the 
ONR  core  program.  The  approach  undertaken  to  solve  the 
issues,  and  the  results  obtained  by  the  investigators  are  sum¬ 
marized.  The  extensive  detailed  information  obtained  from 
this  program  is  omitted  on  purpose.  The  reader  may  refer  to 
the  references  given  at  the  end  of  the  paper  for  further  infor¬ 
mation. 


composites  into  the  supersonic  mixing  layer  and  flame  struc¬ 
ture.  Among  others,  the  program  addressed  the  following 
issues,  which  ate  specific  to  this  paper: 

1.  What  are  the  turbulent  intensities  and  size  scales  in 
compressible  subsonic  and  supersonic  shear  layer  mix¬ 
ing  flows?  Ate  the  structures  large  scale  and  coherent? 

2.  Do  the  stream'vise  narrow  band  velocity  fluctuations, 
momentum  thickness  and  mixedness  of  these  turbulent 
shear  layers  behave  similarly  with  streamwise  distance 

FIgur*  2. 

Temperature  Distribution  in  the  Mixing  Layer  with  and  without 
Energy  Release 


Effect  of  Heat  Release 

To  understand  the  dynamics  and  topology  of  the  coherent 
structures  controlling  the  development  of  reactive  shear  flows, 
in  order  to  develop  both  a  conceptual  model  and  an  analytic 
framework  for  their  description,  Grinstein  perftwmcd  numer¬ 
ical  experiments.  His  approach  is  to  use  direct  and  large-eddy 
numerical  simulation  (LES)  of  spatially  evolving  three-di¬ 
mensional  compressible  shear  flows  with  chemical  reaction 
using  monotone  3D  flux-corrected  transport  (FCT)  algo¬ 
rithms’’^.  The  details  of  the  Reactive  Shear  Flow  Model  and 
the  incorporation  of  chemistry  can  be  found  in  Ref.  3.  The 

Figure  3. 

Mean  Profiles  of  Product  Concentration  and  Temperature  at 
Two  Streamwise  Locations. 
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time-dependent,  compressible,  reactive  flow  conservation 
equations  for  total  mass  and  energy  density,  momentum,  and 
chemical  species  number  densities  are  solved.  The  term  asso¬ 
ciated  with  molecular  diffusion,  viscous  diffusion  and  thermal 
conduction  in  the  energy  equation  are  calculated  explicitly 
using  central-difference  approximations  and  coupled  to  con¬ 
vection  using  timestep- splitting  techniques.  Since  solving  a 
detailed  set  of  chemical  kinetic  rate  equations  in  conjunction 
with  the  unsteady,  three-dimensional  flow  equations  is  beyond 
current  computer  resources,  simplified  global-chemistry  mod¬ 
els  were  used  to  understand  the  effects  of  chemical  energy 
release  on  the  complex  three-dimensional  flow  field. 

Reactive  Shear  Flow  Simulations 

Hussain  and  Hussain  showed  that  coherent  structures  are 
capable  of  enhancing  or  reducing  the  mixing  of  coflowing 
reacting  streams,  and  hence  can  lead  to  methods  of  controlling 
combustion  efficiency  through  control  of  the  initiation,  evolu¬ 
tion,  and  interaction  of  coherent  structures'*.  Chemical  reac¬ 
tions  between  the  coflowing  streams  may  also  affect  the 
structure  of  the  mixing  layer,  either  through  compressibility 
effects  such  as  baroclinic  vorticity  production  and  expansion 
as  energy  is  released,  or  by  changing  the  frequency  content  of 
the  mixing  layers  with  the  introduction  of  characteristic  fre¬ 
quencies  of  the  chemical  processes.  The  combustion  time 
scales  themselves,  in  turn,  may  be  modified  by  the  natural 
acoustic  frequencies  of  the  mixing  layer  or  by  those  imposed 
throu^  excitation.  The  interaction  between  these  effects  can 
have  dramatic  consequences  on  combustion  efficiencies.  Ear¬ 
lier  studies^-^  have  indicted  that  the  energy  released  by  the 
chemical  reaction  has  the  effect  of  reducing  the  entrainment 
in  the  mixing  layer. 

The  results  of  numerical  simulations  of  spatially  evolving, 
chemically  reacting  mixing  layers  ate  presented.  The  system 
studied  consists  of  a  mixing  layer  formed  by  dilute  non- 
premixed  subsonic  streams  of  hydrogen  and  oxygen,  which 
are  allowed  to  react  both  with  and  without  energy  release. 
Instantaneous  results  from  reactive  two-dimensional  mixing 
layer  simulations  are  shown  in  Fig.  1.  Here,  the  product 
concentration  is  plotted  against  the  length  of  the  mixing  layer 
without  energy  release  (Figure  la),  and  with  energy  release 
(Figure  lb),  llie  reaction  considered  is  IHz  +  O2  ■>  2H2O, 
initial  temperature  *  1400  K,  Mach  number  s  0.3,  reactant 
(molar)  concentration  is  10% ,  and  the  velocity  ratio  is  3.  The 
flow  was  organized  by  forcing  the  cross  stream  velocity  at  the 
inflow. 

Spanwise  vortices  rtdl  up  with  the  forcing  frequency  as  a 
result  of  the  nonlinear  evolution  of  the  Kelvin-Helmholtz 
instability,  and  vortex  pairings  are  inhibited.  Thus  the  growth 
of  the  mixing  layer  is  significant  only  in  the  vortex-roll-up 
rqikn.  The  center  of  the  shear  layer  padually  shifts  towards 
the  iqtper,  slower  side  as  the  flow  moves  downstream,  indicat¬ 
ing  tte  irdterently  asymmetric  entrainment  in  the  planar  shear 


nguro  4. 

3-D  Simulation  of  Product  Generation  with  Energy  Release  and 
Spanwise  Excitation. 


:;n.  with  energy  release 


FIgurw  5. 

Enhancement  of  Vorticity  and  Chemical  Reaction  by  Second¬ 
ary  Instabilities 


2-D  Excitation  Only  2-D  +  3-D  Excitation 


VORTICITY  MAGNITUDE 


2-D  Excitation  Only  2-D  +  3-D  Excitation 

INSTANTANEOUS 
PRODUCT  GENERATION 
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flow’’*.  Due  to  the  fast  flow  regimes  considered,  diffusive 
mixing  takes  place  in  the  immediate  neighborhood  of  the 
interface  between  reactants,  and  the  chemical  product  gener¬ 
ation  is  observed  mainly  away  from  the  cores  of  the  structures, 
especially  in  its  outer  edges,  which  are  the  most  chemically 
reactive  regions.  As  can  be  seen  from  Figure  1,  chemical 
energy  release  has  the  effect  of  reducing  the  shear  layer 
growth,  and  the  amount  of  product  formed.  The  product  peak- 
value  is  found  to  be  30%  larger  without  energy  release.  The 
temperature  distribution  with  and  without  eneigy  release  is 
shown  in  Figure  2.  It  is  noted  that  the  chemical  energy  release 
is  accompanied  by  significant  changes  in  the  temperature,  with 
peak-temperatures  at  1.23To  and  1.02To  with  and  without 
energy  release  respectively. 

In  Figure  3,  profiles  of  time-averaged  product  concentra¬ 
tion  and  temperature  are  shown  at  a  streamwise  location  xi 


FIgur*  6. 


Suptrsontc  noixit  Trcvtr«ir>f  probt 


A-A 


Figur*  7. 

Velocity  Profiles  at  Various  Axial  Locations 


Figure  8. 


Experimental  Values  of  Shear  Layer  Spreading  Rate. 
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and  another  streamwise  location  X2  downstream  for  cases  with 
and  without  heat  release.  These  profiles  clearly  show  a  shift 
of  the  flame  towards  its  slower  side,  which  becomes  more 
pronounced  moving  downstream.  Part  of  the  asymmetry  may 
be  attributed  to  the  different  mass  diffusivities  on  each  side  of 
the  shear  layer. 

However,  with  spanwise  excitation,  three-dimensional 
simulations  have  shown  an  increase  in  mixing  even  with  heat 
release.  For  the  same  conditions  used  in  the  jxevious  simula- 
ti(Hi,  the  spatially  evolving  product  concentrations  are  shown 
for  two  consecutive  time  steps  in  Figure  4.  The  three-dimen¬ 
sional  nature  of  the  mixing  layer  and  mixing  enhancement  is 
evident.  Spectral  simulations  were  perfumed  by  Met  Calfe 
and  Hussain  to  understand  the  topology  and  dynamics  of  large 
scale  coherent  structures,  and  their  coupling  with  fine  scale 
mixing.  These  simulations  have  shown  that  secondary  in¬ 
stabilities  enhance  mixing  by  convoluting  and  stretching 
flame  sheet,  and  they  interact  with  “rolls”  to  control  large  scale 
species  transport  and  reaction  rate.  Figure  S  shows  the  vortic- 
ity  magnitude  and  instantaneous  product  generation  with  2-D 
excitation,  and  with  2-D  +  3-D  excitation.  The  three-dimen¬ 
sional  excitation  produces  intense  fine  scale  mixing,  as  shown 
by  the  vorticity  magnitude  and  the  instantaneous  product 
generation’. 

Supersonic  Mixing  Studies 

Mixing  of  Parallel  Supersonic  Streams 

A  study  was  initiated  by  Waltrup  and  Sullins  to  determine 
the  structure  and  growth  of  nonreacting  shear  layers  formed 
when  a  sonic  or  supersonic  two-dimensional  jet  mixes  with  a 
two-dimensional  supersonic  air  stream.  The  primary  airstream 
of  a  typical  Mach  number  of  2  or  3  is  separated  from  a  sonic 
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FIgurm  9. 

Schlieren  Photograph  of  the  Shear  Layer. 
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FlguiwIO. 

Normalized  Growth  Rates  and  Peak  Turbulence  Quantities  as 
a  Function  of  Relative  Mach  Number. 


Mr 


airstream  by  a  sharp  trailing  edge  splitter  plate.  This  is  fol¬ 
lowed  by  a  constant  area  supersonic  mixing  section  (Figure  6). 
Diagnostics  included  Schlieren  photography,  and  pitot  and 
cone  static  pressure  and  velocity  measurements  using  a  one 
component  laser  Doppler  velocimeter. 

If  both  streams  are  incompressible  (M  <  0.3),  the  shear 
layer  is  iidieiently  unstable,  and  laige  scale,  coherent  vortical 
structures  are  formed.  The  growth  is  dependent  on  vortex 
pairing  and  the  ratio  of  the  velocities  of  the  two  streams.  At 
supersonic  q)eeds,  the  shear  layer  tends  to  be  more  stable,  and 
efficient  mixing  is  not  achieved.  The  concept  of  convective 
Mach  number  Me,  established  by  Roshko  and  Papamoschou'” 
applies. 


Tests  were  conducted  with  two-dimensional  shear  layers 
produced  by  both  Mach  2  and  Mach  3  nozzle  streams  mixing 
with  a  Mach  1  nozzle  stream.  Instream  cone  static  and  pitot 
probes  were  used  to  determine  the  velocity  profiles  at  various 
axial  locations  (Figure  7).  The  velocity  profiles  were  used  to 
determine  the  shear  layer  spreading  rate  (Figure  8).  Schlieren 
photographs  were  also  taken  to  determine  the  spreading  rate 
(Figure  9),  and  the  spreading  rate  determined  from  the  velocity 
profiles  compared  well  with  these  photographs. 

In  a  series  of  experiments,  Krier  and  his  associates  inves¬ 
tigated  the  fundamental  fluid  dynamic  and  combustion  mech¬ 
anisms  occurring  in  the  entrainment,  mixing,  and  combustion 
processes  in  the  free  shear  layer  formed  between  two  high 
speed  (supersonic/supersonic  or  supersonic/subsonic) 
streams". 

Mean  and  Turbulent  Velocity  Measurements:  Velocity 
measurements  have  been  obtained  using  a  two-component 
laser  Doppler  velocimeter  system  in  a  two  stream,  supersonic 
mixing  layer  facility.  A  variety  of  cases  has  been  examined 
with  free  stream  velocity  ratios  in  the  range  from  0.16  toO.78. 
freestream  density  ratios  ranging  from  0.S7  to  1.SS,  and  rela¬ 
tive  Mach  numb^  from  0.40  to  1.97.  (The  relative  Mach 
number  is  used  here  as  a  measure  of  compressibility  effects 
and  is  defined  as  the  freestream  velocity  difference  divided  by 
the  mean  speed  of  sound,  AftS  A  U/ a).  The  range  of  relative 
Mach  numbers  examined  spans  the  range  of  significant  com¬ 
pressibility  effects  as  determined  by  the  reduction  in  normal¬ 
ized  shear  layer  growth. 

Measurements  from  the  developing  regions  of  these  shear 
layers  suggest  that  a  local  Reynolds  number  based  on 
freestream  velocity  difference  and  mixing  layer  thickness  of 
Rcfd  s  A  UB/  V  =  1  X  10*  is  required  for  full  development  of 
the  mean  turbulent  velocity  fields.  To  determine  the  growth 
rates,  least-squares  linear  regressions  were  performed  on  the 
mixing  layer  thickness  data  from  the  fully  developed  region. 
The  resulting  growth  rates  have  been  normalized  by  the 
growth  rate  of  corresponding  incompressible  mixing  layer  at 


Flgurall. 

2-D  Scalar  Covariance  Field  for  Mr  =  0.63.  (a)  Trivisverse 
Image,  (b)  Oblique  Image 
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the  same  freestieam  velocity  and  density  ratios,  and  these 
normalized  growth  rates  are  plotted  in  Figure  10  with  respect 
to  the  relative  Mach  number.  The  results  clearly  show  a 
reduction  in  the  normalized  growth  rate  due  to  compressibility, 
which  is  consistent  with  those  observed  by  Roshko  and  others. 

Variations  in  the  peak  values  of  the  fully  developed  high¬ 
speed  mixing  layer  turbulence  quantities  from  their  incom¬ 
pressible  counterparts  reflect  the  effects  of  compressibility  on 
turbulence.  The  peak  values  of  the  streamwise  turbulence 
intensity,  transverse  turbulence  intensity,  and  the  normalized 
kinematic  Reynolds  shear  stress,  -<u'  v'>/(  AU)^  for  the  com¬ 
pressible  cases  have  been  normalized  by  the  respective  incom¬ 
pressible  values  of  0.18,  0.14,  and  0.013  and  are  plotted  in 
Figure  10.  From  order-of-magnitude  estimates  using  boundary 
layer  approximations  to  the  governing  equations,  it  can  be 
shown  that  for  a  compressible  mixing  layer  the  normalized 
kinematic  Reynolds  shear  stress  should  scale  with  the  normal- 


Flgural  2. 

2-D  Scalar  Covariance  Field  for  Mr  =  1.49.  (a)  Transverse 
Image,  (b)  Oblique  Image. 
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Comparison  of  Product  Generated  with  and  without  3-D  Exci¬ 
tation. 
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ized  growth  rate.  The  comparable  reduction  in  these  two 
quantities  with  increasing  relative  Mach  number  is  seen  in 
Figure  10. 

Large  Scale  Structures:  The  large  scale  structures  of  com¬ 
pressible  mixing  layers  have  been  studied  by  Krier’s  group 
using  planar  Mie  scattering  imaging  and  statistical  analysis  of 
two-dimensional  covariance  fields’^.  Three  flow  conditions 
were  investigated,  using  both  traditional  transverse  imaging 
as  well  as  oblique  imaging.  Two  types  of  seeding  methods 
were  employed,  passive  scaler  transport  of  sub-micron  ethanol 
droplets  from  one  freestream  to  the  other  as  well  as  the 
formation  of  ethanol  droplets  due  to  molecular  mixing  of  air 
from  one  stream,  laden  with  ethanol  vapor,  with  cold  air  from 
the  other  stream.  The  relative  Mach  numbers  of  the  three  flow 
conditions  were  M,  =  0.63,  1.24  and  1.49,  while  the  local 
Reynolds  numbers  for  these  cases  were  Re  =  A  Ub/  v  =  3.0 
X 10*  and  6.3  x  lO*  respectively.  The  level  of  three-dimension¬ 
ality  in  a  compressible  mixing  layer  was  found  to  be  largely  a 
function  of  relative  Mach  number  with  the  Reynolds  number 
^parently  having  little  effect  other  than  to  introduce  a  wider 
range  of  turbulent  structure  scales. 

The  instantaneous  images  provide  some  striking  exam¬ 
ples  of  the  various  flow  structures  within  compressible  mixing 
layers.  Organized  large  scale  structures  were  present  in  all 
cases;  however  they  appeared  most  frequently  in  the  low 
relative  Mach  number  flow.  As  the  relative  Mach  number  is 
increased,  the  transverse  images  suggest  the  structures  become 
less  organized,  with  an  increased  presence  of  randomly  ori¬ 
ented,  smaller  scale.  Structures  indicative  of  streamwise 
counter-rotating  vorticity,  and  a  measure  of  the  level  of  three-di¬ 
mensionality  of  the  mixing  layer,  were  more  often  observed  at 
higher  relative  Mach  number  cases.  The  sizes  of  these  structures 
also  tended  to  increase  with  the  frequency  of  occurrence. 


Figur«14. 

Typical  Convoluted  Splitter  Surface  Showing  Periodic  Stream- 
wise  Vortex  Array. 
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The  spatial  covariance  flelds  of  the  ensembles  of  2S6 
instantaneous  images  for  each  case  were  computed  to  deter¬ 
mine.  on  a  sound  statistical  basis,  the  characteristics  of  the 
large  scale  structures.  Figure  11  and  12  show  the  covariance 
fields  for  both  the  transverse  and  oblique  scalar  transport 
images  of  the  low  and  high  relative  Mach  number  cases, 
respectively.  In  Figure  11  and  12,  the  central  peak  is  100%, 
and  each  successive  contour  is  20%  lower.  These  results 
indicate  that,  as  the  relative  Mach  number  is  increased,  the 
transverse  scalar  transport  structures  tend  to  increase  in 
dimensionless  size  ratio  and  eccentricity  while  the  angular 
orientation  decreases  slightly.  These  results  agree  extremely 
well  with  recent  time-evolving  numerical  simulations’^  of 
comiHessible  mixing  layers  which  also  show  more  flattened, 
elliptical  structures  at  large  relative  Mach  number.  Oblique 
scalar  transport  correlations  indicate  structures  of  generally 
smaller  size  ratio  and  significantly  reduced  eccentricity  than 
their  transversely  imaged  counterparts,  yet  the  effect  of  in¬ 
creased  relative  Mach  number  indicates  only  a  slight  increase 
in  dimensionless  structure  size  with  no  apparent  trend  in  the 
eccentricity.  The  product  formation  studies  exhibited  struc¬ 
tures,  in  both  the  transverse  and  oblique  image  planes,  which 
were  smaller  in  dimensionless  size  and  less  elliptical  in  shape 
than  the  top  stream  scalar  transport  structures.  Significantly 
different  characteristics  between  transverse  and  oblique  real¬ 
izations  of  the  scalar  transport  structure  implies  the  presence 
of  an  organized,  large  scale  mixing  layer  structure,  which  is 
not  entirely  due  to  increased  three-dimensionality  of  the  mix¬ 
ing  layer  associated  with  increased  relative  Mach  number. 


FiguiwIS. 

Comparison  of  Flame  Spreading  Achieved  with  a  Flat  Splitter 
and  a  Convoluted  Splitter. 


Figur*16. 

Vortex  Ring  Penetration  with  Transverse  Pulse  Injection 


Pitot  Tube  Growth  with  Planar  Shock  Impingement  on  bound¬ 
ary  Layer. 
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Mixing  and  Combustion 
Enhancement 

Subsonic  Flows 

As  mentioned  in  section  II,  three-dimensional  excitation 
has  been  found  to  increase  mixing  and  combustion.  The  en¬ 
hancement  in  vorticity  and  product  generation  by  excitation, 
shown  in  Figure  5,  is  Airther  quantiHed  in  Figure  13.  Here  the 
ratio  of  product  generation  with  3-D  excitation  to  that  with  2-D 
excitation  is  plotted  against  time.  As  can  be  seen,  over  30% 
increase  in  product  genoation  is  obtained  by  3-D  excitation. 

Streamwise  Vorticity  Stirring:  Experimental  work  by 
McVey  focussed  on  the  effect  of  streamwise  vorticity  on  the 
spreading  of  flames  in  high  speed  flows.  Using  a  convoluted 
splitter  surface  as  shown  in  Figure  14 ,  and  by  the  proper  choice 
of  the  amplitude,  pitch  and  profile  of  the  leered  surface,  it  was 
possible  to  generate  arrays  of  very-large-scale,  intense  vorti¬ 
ces.  Flow  visualizations  indicated  that  the  influence  of  these 
vortices  is  first,  to  cause  a  large  scale  exchange  of  fluid 
(stirring)  followed  by  a  rapid  breakdown  of  the  organized 
strtictures  into  random  small  scale  turbulence,  which  is  essen¬ 
tial  to  high  speed  combustion 

Non-premixed  flame  with  air  as  the  oxidant  and  a  mixture 
of  nitrogen  and  vaporized  triethylborane  as  fuel  was  used. 
Since  the  fuel  is  pyrophoric,  no  ignition  source  was  required. 
Flame  photographs  with  a  flat  splitter  and  a  lobed  splitter  are 
shown  in  figure  IS.  A  dramatic  increase  in  the  rate  of  flame 
spreading  is  evident,  which  is  attributed  to  the  streamwise 
vorticity  stirring.  Further  work  to  understand  the  details  of  the 
mechanisms  involved  is  presently  sponsored  by  the  Army 
Research  Office. 


FigurelS. 

Vortex  Generators 


(a)  Circular  cylinders 


Fig.  19  Vortex  Genmton 


Figura20. 


FlguMlS. 


Pitot  Thickness  Growth  with  Cylindrical  Vortex  Generators. 
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Transverse  Pulsed  Fuel  Injection:  Vakili  and  Wu  are  sttidy  ing 
transverse  pulse  injection  of  fuel  into  the  main  oxidant  stream,  as 
a  means  of  enhancing  mixing  and  combustion.  Experiments 
conducted  in  subsonic  flows  indicated  that  significant  penetration 
of  the  vortex  ring  can  be  obtained  by  transverse  pulse  injection 
(Figure  16).  In  particular,  sharp  edge  orifices  have  been  found  to 
be  very  effective.  Further  studies  will  be  undertaken  to  determine 
the  effectiveness  in  super  onic  mixing  enhancement,  and  conse¬ 
quently  supersonic  combustion'^. 

Supersonic  Flows 

Compressible  shear  layers  are  very  stable,  and  the  spread¬ 
ing  rate  decreases  as  the  convective  Mach  number  increases. 
In  order  to  achieve  complete  combustion  within  the  short 
residence  times  available  in  supersonic  combustors,  it  is  nec¬ 
essary  to  enhance  mixing.  Energy  release  could  further  com¬ 
plicate  and  decrease  mixing.  Experimental  studies  have  been 
undertaken  by  Dolling,  and  by  Schadow  to  enhance  supersonic 
mixing,  and  mixing  and  combustion  respectively,  by  passive 
or  active  means. 

Dolling’s  research  focused  on  the  growth  rate  and  struc¬ 
ture  of  high  Reynold’s  number  turbulent  shear  layer  formed 
by  Mach  5  and  Mach  3  airstreams'*.  The  first  phase  of  the 
study  was  to  determine  whether  the  growth  rate  of  the  shear 
layer  could  be  increased  by  changing  the  initial  conditions  at 
the  shear  layer  origin.  The  initial  conditions  were  changed 
either  by  a  planar  3-D  shock  impingement  on  the  turbulent 
boundary  layer  upstream  of  the  origin,  or  by  using  cylindrical 
or  wedge  shaped  vortex  generators  in  the  boundary  layer 
upstream  of  the  origin. 

Shock  Impingement:  Both  planar  and  3-D  shock  waves 
were  generated  and  arranged  to  impinge  on  either  the  bound¬ 
ary  layer  upstream  of  the  shear  layer  origin  or  at  the  shear  layer 
itself.  Seven  cases  were  studied.  Typical  results  will  be  pre¬ 
sented  here. 

For  planar  shock  impingement  on  boundary  layer,  three 
different  configurations  were  used.  In  Figure  17,  the  growth 


rate  of  the  shear  layers  for  these  cases  are  shown.  The  growth 
rate  of  the  undisturbed  shear  layer  is  also  shown,  by  hatch  line, 
for  comparison.  A  maximum  of  30%  increase  in  shear  layer 
growth  rate  was  obtained.  However,  within  a  few  boundary 
layer  thicknesses  downstream,  the  growth  rates  are  not  much 
different  from  the  undisturbed  growth  rate. 

Three-dimensional  shock  impingement  on  the  shear  layer 
generated  only  slight  three-dimensionality  in  the  layer  im¬ 
mediately  downstream  of  the  impingement  location  due  to  the 
pressure  gradient  caused  by  the  sine  wave  shape  of  the  incident 
^ock.  However,  the  flow  became  mwe  two-dimensional  fur¬ 
ther  downstream,  and  the  thickness  of  the  shear  lay^  remained 
about  the  same.  The  growth  rates  across  the  span  are  summa¬ 
rized  in  Figure  18.  As  can  be  seen,  th^e  are  no  observable 
differences  between  the  growth  rate  and  thickness  for  this  case 
with  the  undisturbed  case. 

Vortex  Generation:  Two  types  of  vortex  generators  were 
used.  The  first  consisted  of  a  row  of  small  cylinders  screwed 
into  the  test  surface,  and  the  second  consisted  of  corrugated 
plates  or  wedges  (Figure  19).  The  results  show  that  the  shear 
layer  growth  rate  can  be  increased  by  about  30%  using  cylin¬ 
drical  vortex  generator  arrangements  (Figure  20),  whereas 
with  wedge-shaped  arrangements,  the  growth  rate  is  unaf¬ 
fected  or  even  reduced. 

Non  axisymmetric  Nozzle  Flows:  Schadow  and  his  asso¬ 
ciates  have  demonstrated  the  superiority  of  non  axisymmetric 
nozzles,  with  respect  to  circular  nozzles,  in  terms  of  fine  scale 
and  large  scale  mixing  in  free-jet  and  shrouded-jet  tests  with 
and  without  reaction'^’’*.  They  have  extended  this  work  to 
supersonic  jets.  The  experiments  were  performed  in  a  coaxial 
supersonic  jet  facility,  which  can  be  used  for  nonreacting 
mixing  studies  and  for  supersonic  fuel-rich  plume  combustion 
studies  (Figure  21). 

For  the  mixing  studies,  the  inner  (center)  jet  was  kqrt  at 
Ml  =  3.  Helium,  argon,  nitrogen,  and  heated  and  unhealed  air 
were  used  to  change  the  density  of  the  inner  jet.  The  velocity 

Figi<ra22. 


FlgufwZI. 


Normalized  Spreading  Rate  of  Circular  and  Noncircular  Jets 
as  a  Function  of  Convective  Mach  Number. 
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Pigura23. 

Radiance  Contours  of  Circular  and  Multi-step  Nozzles. 


Figur«24. 

Relative  Total  Radiance  for  Circular  and  Multi-step  Nozzles. 


RELATIVE 

RADIANCE 


CIRCULAR  MULTI-STEP 


of  the  outer  jet  was  varied  from  zero  (no  coaxial  flow)  to  fully 
expanded  condition  (M2  =  1.8)  for  ambient  and  47SK  air 
temperatures.  Six  test  conditions  were  selected  to  obtain  a 
convective  Mach  Numbo-  range  of  0.5  Me  2.2.  Mixing  be¬ 
tween  the  center  jet  and  outer  air  was  determined  with  a 
total-pressure  probe,  and  gas  sampling  probe  connected  to  an 
analyzer. 

In  Figure  22,  the  variation  of  the  normalized  spreading 
rates  of  the  circular  jet  and  the  rectangular  jets  (two  conHgu- 
rations)  are  shown  against  the  convective  t^h  number.  The 
behavior  observed  for  the  circular  coaxial  jet  is  similar  to  that 
of  the  planar  shear  layer,  with  the  normalized  qveading  rate 
decreasing  with  increasing  convective  Mach  number.  How¬ 
ever,  the  qneading  rates  of  the  rectangular  jets  are  higher  than 
the  circidar  jet,  both  in  the  major  and  minor  axis  planes,  in 
almost  the  entire  convective  MKh  number  range.  The  spread¬ 


ing  rates  are  relatively  higher  for  Me  0.4,  because  the  rectan¬ 
gular  jets  have  a  higher  incompressible  spreading  rate  than  the 
circular  jet 

The  supersonic  fuel-rich  plume  combustion  experiments 
were  also  conducted  in  the  same  coaxial  jet  facility.  Here,  the 
outer  flow  at  ambient  temperature  was  fixed  at  Me  =  1.3  and 
the  inner  fuel-rich  flow  at  Mi  =  2.  The  fuel-rich  plume  was 
{X'oduced  in  a  gas  generator  by  burning  ethylene,  air,  and 
oxygen  at  an  equivalence  ratio  of  2,  and  a  combustion  temper¬ 
ature  of  Ti  =  2300K.  Combustion  characteristics  are  compared 
based  on  the  infrared  images,  which  indicate  the  spreading  rate 
and  combustion  intensity.  As  a  typical  example,  in  Figure  23, 
the  radiance  map  of  the  underexpanded  coaxial  jets  from  a 
circular  and  multi-step  nozzle  are  compared  ’  ^  The  RMS  level 
of  the  multi-step  nozzle  is  higher  than  the  circular  jet  In  Figure 
24,  the  circular  and  multi-step  nozzles  are  compared  for  three 
{nessure  ratios  based  on  the  relative  integrated  radiance  from 
the  nozzle  lip  to  x/De  =  16.  The  radiance  is  normalized  by  the 
value  for  the  circular  jet  at  the  design  pressure  ratio.  These 
figures  indicate  the  superior  performance  of  the  multi-step 
nozzle,  and  the  increased  heat  release  associated  with  en¬ 
hanced  fine-scale  mixing. 


Summary 

A  special  focus  program  was  sponsexed  by  ONR  to  un¬ 
derstand  the  mechanisms  involved  and  the  control  methodol¬ 
ogy  to  enhance  mixing  and  combustion  is  subsonic  and 
supersonic  reacting  and  nonreacting  flows.  This  multi-disci¬ 
plinary,  numerical  and  experimental  research  by  the  academia. 
Navy  and  industry  laboratories  produced  the  following  signif¬ 
icant  results: 

1.  Energy  release  has  a  tendency  to  reduce  mixing  and 
further  product  generation. 

2.  Three-dimensional  excitation  increases  product  genera¬ 
tion  (with  energy  release). 

3.  Secondary  instabilities  are  very  important  in  both  the 
core  and  braid  regions  of  the  flow,  and  they  interact  with 
the  rolls  to  control  large  scale  species  transport,  and  thus 
the  reaction  rate. 

4.  Streamwise  vorticity  stirring  by  means  of  convoluted 
splitter  plates  has  shown  a  significant  improvement  in 
flame  spreading. 

5.  Transverse  pulsed  injection  of  the  fuel  allows  deep 
penetration  of  the  vortex  rings  into  the  mainstream,  and 
enhances  mixing. 

6.  Supersonic  mixing  layers  are  very  stable,  and  the  growth 
of  the  mixing  layer  decreases  as  the  convective  Mach 
number  increases. 

7.  Shock  impingement  in  the  boundary  layer  in  supersonic 
flow  produces  a  small  increase  in  the  shear  layer  growth , 
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but  within  a  few  boundary  layer  thicknesses  down¬ 
stream,  it  becomes  similar  to  the  undisturbed  growth. 

8.  Supersonic  shear  layer  growth  is  increased  by  about 
30%  by  vortex  generation,  using  cylindrical  vortex  gen¬ 
erators. 

9.  Non  axisymmetric  nozzles  have  increased  mixing  rate 
as  compared  to  circular  nozzles.  In  coaxial  supersonic 
flows,  they  have  produced  increased  flne  scale  mixing, 
and  enhanced  combustion. 

Further  research  and  develqrment  is  needed  to  control 
and  enhance  supmonic  mixing  and  combustion.  Active,  pas¬ 
sive,  and  hybrid  control  methodologies  may  be  required  fcx- 
specific  application.  OMR  is  currently  pursuing  another  spe¬ 
cial  focus  program  on  active  control  of  complex  flows. 
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Abstract 

A  fcnmal  approach  for  evaluating  the  reliability  of  com¬ 
puter  software  is  through  probabilistic  models  and  statistical 
methods.  This  paper  is  an  expository  overview  of  the  literature 
on  the  former.  The  various  probability  models  for  software 
failure  can  be  classiHed  into  two  groups;  the  merits  of  these 
groups  are  discussed  and  an  example  of  their  use  in  decision 
problems  is  given  in  some  detail.  The  direction  of  current  and 
future  research  is  contemplated.  This  research  has  been  sup¬ 
ported  by  the  Office  of  Naval  Research,  the  Army  Research 
Office,  and  the  Air  Force  Office  of  Scientific  Research. 

1.  Introduction 

Having  been  developed  over  the  last  20  years,  software 
reliability  is  a  relatively  new  area  of  research  for  the  statistics 
and  the  computer  science  communities.  It  arose  because  of 
interest  in  trying  to  predict  the  reliability  of  software,  particu¬ 
larly  when  its  failure  could  be  catastrophic.  Obviously  the 
software  that  controls  an  aircraft  carrier,  a  nuclear  power 
station,  a  submarine  or  a  life-suj^rt  machine  needs  to  be  very 
reliable,  and  statistical  techniques  will  aid  the  computer  scien¬ 
tist  in  deciding  if  such  software  has  sufficient  reliability.  The 
subject  is  also  of  commercial  importance,  as  for  example  when 
decisions  have  to  be  made  concerning  the  release  of  software 
into  the  marketplace. 


All  software  is  subject  to  failure,  due  to  the  inevitable 
presence  of  errors  (or  bugs)  in  the  code,  so  the  first  aim  of  the 
subject  has  been  to  develop  models  that  describe  software 
failure.  There  are  various  methods  of  specifying  such  failure 
models,  and  Section  2  discusses  these  in  some  detail.  It  is  fair 
to  say  that  this  model  derivation  has  been  the  focus  of  research 
so  far.  Once  a  failure  model  has  been  specified  then  it  can  be 
applied  to  problems  such  as  the  optimal  time  to  debug  software 
or  deciding  whether  software  is  ready  for  release.  These  appli¬ 
cations  have  received  less  attention  in  the  literature  but  are 
becoming  more  prevalent  We  will  mention  here  that  there  is 
another  approach  to  software  reliability  that  differs  consider¬ 
ably  from  the  statistical  ideas  presented  here.  This  approach 
attempts  to  prove  the  reliability,  or  correctness,  of  software  by 
formal  means  of  proof,  just  as  one  would  prove  a  mathematical 
theorem.  This  is  an  exercise  in  logic,  albeit  a  rather  complex 
one.  It  works  well  on  small  programs,  for  example  on  a 
program  that  computes  the  factorial  function,  but  becomes  a 
forbidding  task  for  even  moderately  complex  pieces  of  code. 
Nevertheless,  the  idea  that  software  can  be  proved  correct  is 
appealing.  The  approach  is  not  discussed  further. 

This  paper  is  divided  into  5  further  sections.  Section  2 
categorizes  the  different  strategies  that  have  been  used  to 
model  software  failure.  Section  3  reviews  the  historical  devel¬ 
opment  of  the  subject  by  describing  some  of  the  more  com¬ 
monly  used  models,  and  Section  4  shows  that  many  of  these 
models  can  be  unified  if  one  adopts  a  Bayesian  position. 
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FIgur*  1. 

Categorization  of  Software  Reliability  Models. 


Section  S  looks  at  applications  of  the  material  developed  in 
Section  3,  and  Section  6  concludes  with  a  look  at  the  current 
and  future  direction  of  the  subject  We  assume  that  the  reader 
has  some  familiarity  with  some  basic  reliability  and  probabil¬ 
ity  concepts;  in  particular  it  is  important  that  he  or  she  has 
knowledge  of  some  common  probability  distributions,  statis¬ 
tical  inference  and  decision  making,  Poisson  processes  and  the 
concept  of  a  failure  rate. 


2.  Model  Categorization 

All  statistical  software  reliability  models  are  probabilistic 
in  nature.  They  attempt  to  specify  the  probability  of  software 
failure  in  some  manner.  In  looking  thrwgh  the  literature,  one 
observes  that  the  models  developed  so  far  can  be  broadly 
classified  into  two  categories. 

Type  I:  Those  which  propose  a  probability  model  for 
limes  between  successive  failure  of  the  software,  and 

Type  II:  Those  which  prcqMse  a  probability  model  for 
the  number  up  to  a  certain  time.  Time  is  often  taken 

to  be  CPU  time,  or  the  amount  of  time  that  the  software  is 
actually  ranning,  as  opposed  to  real  time.  In  theory,  specifica¬ 
tion  via  one  of  these  two  methods  enables  one  to  specify  the 
other.  So  a  model  that  specifies  time  between  failure  will  also 
be  able  to  tell  you  the  number  of  failures  in  a  given  time,  and 
vice  versa,  hi  practice,  this  may  not  be  straightforward. 

The  fiist  of  these  categories,  modeling  time  between 
failure,  is  roost  commonly  accomplished  via  a  specification  of 
the  failure  rate  of  the  software  as  it  is  running.  When  this  is 
the  case  the  model  is  to  be  of  Type  I- 1 .  The  failure  rate  for  the 
i-th  time  between  failure  is  given,  for  i^l,  2,  3,...  and  a 
probability  modd  results.  One  distinctive  feature  of  software 


is  that  its  failure  rate  may  decrease  with  time,  as  more  bugs 
are  discovered  and  corrected.  This  contrasts  with  most  me¬ 
chanical  systems  which  will  age  over  time  and  so  have  an 
increasing  failure  rate.  An  attempt  to  debug  software  may 
introduce  more  bugs  into  it,  thus  tending  to  increase  the  failure 
rate,  so  the  decreasing  failure  rate  assumption  is  somewhat 
idealized.  However,  most  of  the  models  of  this  type  that  are 
reviewed  here  have  a  decreasing  failure  rate. 

Another  way  to  model  time  between  failure  is  to  define  a 
stochastic  relationship  between  successive  failure  times.  Mod¬ 
els  that  are  specified  by  this  method  are  known  as  Type  1-2, 
and  have  the  advantage  over  Type  I-l  in  that  they  model  the 
times  between  failure  directly,  and  not  via  the  abstract  concept 

of  a  failure  rate.  For  example,  let  Tj ,  Tj . T, ...  be  the  length 

of  time  between  successive  failure  of  the  software.  As  a  simple 
case,  one  could  declare  that  T^j  =  pT  +  e^,  where  p  >  0  is  a 
constant  and  e.  is  an  error  term  (typic^ly  some  random  vari¬ 
able  with  mean  0).  Then  p  <  1  would  indicate  decreasing  time 
between  failure  (software  reliability  expected  to  become 
worse),  p  =  1  would  indicate  no  expected  change  in  software 
reliability  whilst  p  >  1  indicates  increasing  times  between 
failure  (software  reliability  expected  to  improve).  Those  fa¬ 
miliar  with  time  series  will  recognize  the  relationship  in  this 
example  as  an  auto-regressive  process  of  order  1;  in  general, 
one  would  say  T^j  =  f(Tj,  Tj, ...,  T.)  +  for  some  function  f. 

The  second  of  these  categories,  modeling  the  number  of 
failures,  uses  a  point  process  to  count  the  failures.  Let  M(t)  be 
the  number  of  failures  in  the  software  that  are  observed  during 
time  [0,t).  M(t)  is  modeled  by  a  Poisson  process,  which  is  a 
stochastic  process  with  the  following  properties. 

i)  M(0)  =  0  and  if  s  <t  then  M(s)  S  M(t).  M(t)'  takes 
values  in  {0, 1,2,...) 
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ii)  The  number  of  failures  that  occur  in  disjoint  time  inter¬ 
vals  are  independent.  So,  for  example,  the  number  of 
failures  in  the  first  S  hours  of  use  has  no  effect  on  the 
number  of  failures  in  the  next  5  hours. 

iii)  The  number  of  failures  to  time  t  is  a  Poisson  random 
variable  with  mean  |j.(t),  for  some  non-decreasing  func¬ 
tion  p.(t);  that  is  to  say: 

P(M(t)  =  n)  =  ^^^^e-'‘W  n  =  0,1,2,... 
n! 

The  different  models  of  this  type  have  a  different  function 
p(t),  which  is  called  the  mean  value  function.  The  mean  number 
of  failures  at  time  t  is  indeed  p(t),  as  is  the  variance.  The  Poisson 
[Htxess  is  chosen  because  in  many  ways  it  is  the  simplest  point 
{Hocess  yet  it  is  flexible  and  has  many  useful  properties  that  can 
be  exploited.  This  second  approach  has  become  increasingly 

Figur*  2(1). 

The  failure  rate  of  the  model  of  Jelinski  ar)d  Moranda 


r(t) 


Flgur*  2(11). 

The  failure  rate  of  the  model  of  Littlewood  and  Verall. 


Hi) 


popular  in  recent  years.  M(t)  can  also  be  specified  by  its 
intensity  function  X(t),  which  is  the  derivative  of  p(t)  with 
respect  to  t;  either  of  these  functions  completely  qiecify  a  partic¬ 
ular  Poisson  process.  One  disadvantage  of  this  ai^iroach  is  that  it 
implies  that  there  are  conceptually  an  infinite  number  of  bugs  in 
the  program,  which  is  obviously  impossible  for  code  of  a  finite 
length.  Another  disadvantage  is  more  subtle;  the  model  implies  a 
positive  correlation  between  the  number  of  failures  in  adjoining 
time  intervals,  a  situation  which  is  not  true  since  again  the  total 
number  of  bugs  has  to  be  finite. 

Figure  1  is  a  flow-chart  showing  the  above  categorization 
of  the  statistical  models. 

3.  Review  of  Some  Software 
Reliability  Models 

This  section  introduces  some  of  the  more  well  known 
probability  models  for  software  reliability.  There  are  examples 
of  models  from  each  of  the  two  main  categories  that  were 
discussed  in  the  previous  section.  Since  the  main  purpose  of 
the  review  is  to  describe  the  ideas  and  assumptions  behind  the 
models,  technical  details  will  be  kept  to  a  minimum  in  most 
cases.  Those  interested  in  the  details  of  a  particular  model  are 
advised  to  reference  the  papers  where  they  were  originally 
presented. 

Some  common  notation  will  be  assumed  throughout  this 
section  and  is  given  below; 

i)  Ti=  i-th  time  between  failure  of  the  software  [i.e.  time 
between  (i-l)th  and  i-th  failure]. 

ii)  'Ti®=  failure  rate  for  Ti,  the  i-th  time  between  failure,  at 
timet 

iii)  M(t)=  number  of  failures  of  the  software  in  the  time 
interval  [0,  t)  (a  Poisson  process). 

Figur*  2(iii) 

The  failure  rate  of  the  model  of  Moranda. 

Ml)  II 
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iv)  X(t)=  intensity  function  of  M(t). 

v)  ti(t)=  expected  number  of  failures  of  software  in  time 

[0.0. 

=  X  (s)  ds,  since  M(t)  is  a  Poisson  process. 

10  models  are  presented.  Model  numbers  1  to  7  are  of 
Type  I-l,  models  8  and  9  are  of  Type  II  and  model  10  is  of  type 
1-2.  A  common  problem  to  all  the  models  is  the  lack  of  data  on 
which  to  test  their  validity;  data  on  software  reliability  is 
commercially  sensitive  and  so  statisticians  in  academia  have 
very  little  information  on  how  software  in  the  marketplace 
actually  performs.  For  this  reason  it  is  important  that  the 
assumptions  made  in  deriving  these  models  are  carefully 
thought  about 

The  Model  of  Jellnski  &  Moranda 
(1972). 

This  was  the  very  first  software  reliability  model  that  was 
pressed,  and  has  formed  the  basis  for  manv  models  devel¬ 
oped  after.  It  is  a  Type  I-l  model;  it  models  times  between 
failure  by  considering  their  failure  rates.  Jelinski  and  Moranda 
reasoned  as  follows.  Suppose  that  the  total  number  of  bugs  in 
the  program  is  N,  and  suppose  that  each  time  the  software  fails, 
one  bug  is  corrected.  The  failure  rate  of  the  i-th  time  between 
failure,  Tj,  is  then  assumed  a  constant  prc^rtional  to  N-i+1, 
which  is  the  number  of  bugs  remaining  in  the  program.  In 
other  words 

r.j  (f  I N.  A)  =  A  (N-i  +1),  i  =  1,2,3,...,/  >  0. 
for  some  constant  A  . 


Figura  2(lv). 

The  failure  rate  of  the  model  of  Schick  and  Wolverton. 
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There  are  some  criticisms  that  one  could  make  of  the 
model.  It  assumes  that  each  error  contributes  the  same  amount 
to  the  failure  rate,  whereas  in  reality  different  bugs  will  have 
different  effects.  It  also  assumes  that  every  time  a  fix  is  made, 
no  new  bugs  are  introduced;  note  [see  Figure  2(i)]  that  the 
successive  failure  rates  are  indeed  decreasing.  A  model  like 
this  is  sometimes  referred  to  as  a  "de-eutrophication  model", 
because  the  process  of  removing  bugs  from  software  is  akin 
to  the  removal  of  pollutants  in  rivers  and  lakes. 

Bayesian  Reliability  Growth  Model 
(Littlewood  &  Verall  (1973)). 

Like  the  Jelinski  &  Moranda  model,  the  model  proposed 
by  Littlewood  and  Verall  looked  at  times  between  failure  of 
the  software.  However,  they  did  not  develop  the  model  by 


Figure  2(v). 

The  intensity  function  for  the  model  of  Goel  and  Okumoto. 


Figure  2(vi). 

The  intensity  function  for  the  model  of  Musa  and  Okumoto. 
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characterizing  the  failure  rate;  rather  they  stated  that  the  model 
should  not  be  based  on  fault  content  (as  Jelinski  &  Moranda 
had  assumed)  and  then  declared  that  T.  has  an  exponential 
distribution  with  scale  A.,  and  that  A.  itself  has  a  gamma 
distribution  with  shape  a  and  scale  \|/(i),  for  some  function  \|/. 
Despite  this  it  is  still  considered  to  be  a  IVpe  I-l  model. 

Specifically: 


(t  I  A.)  =  A.  e-V  tSO 


n  .  =  X>0 

i  f  (a) 

y(i)  was  supposed  to  describe  the  quality  of  the  program¬ 
mer  and  the  programming  task.  As  an  example,  they  chose 
y(i)=  P|i.  One  can  show  that  this  makes  the  failure  rate 
of  each  T.  decreasing  in  t  and  that  each  time  a  bug  is  discovered 
and  fixed  there  is  a  downward  jump  in  the  successive  failure 
rates:  see  Figure  2(ii).  In  fact 


If  P,  >I  then  the  jumps  in  the  failure  rate  decrease  in  i,  if 
P,  <1  they  increase  whilst  if  P,  =  1  they  remain  a  constant.  So 
if  3,  differs  from  1  then  the  fixing  of  each  bug  is  making  a 
different  contribution  to  the  reduction  in  the  failure  rate  of  the 
software,  an  apparent  advantage  over  the  model  by  Jelinski  & 
Moranda.  This  model  has  received  quite  a  lot  of  attention  and 
has  been  the  subject  of  various  modifications;  see  models  6 
and  7  later  in  this  section. 


The  De-Eutrophication  model  of 
Moranda  (1875). 

Another  (de-eutrophication)  model  of  Moranda  (1975) 
attempted  to  answer  some  of  the  criticisms  of  the  Jelinski  & 
Moranda  model,  in  particular  the  criticism  concerning  the 
equal  effect  that  each  bug  in  the  code  has  on  the  failure  rate. 
He  hypothesized  that  the  fixing  of  bugs  that  cause  early 
failures  in  the  system  reduces  the  failure  rate  mote  than  the 
fixing  of  bugs  that  occur  later,  because  these  early  bugs  are 
more  likely  to  be  the  bigger  ones.  With  this  in  mind,  he 
proposed  that  the  failure  rate  should  remain  constant  for  each 
T.,  but  that  it  should  be  made  to  (tectease  geometrically  in  i 
after  each  failure  i.e.  for  constants  D  and  k. 

r.^  (tlDJc)  =  Dk^'  t^O,D>0andO<k<l. 


Imperfect  Debugging  Model 
(Goel  &  Okumoto  (1978)). 

This  model  is  another  generalization  of  the  Jelinski  & 
Moranda  model  which  attempts  to  address  the  criticism  that  a 
perfect  fix  of  a  bug  does  not  always  occur.  Goel  &  Okumoto 's 
Imperfect  Debugging  Model  is  like  the  Jelinski  &  Moranda 
model,  but  assumes  that  there  is  a  probability  p,  0  <  p  <  1 ,  of 
Hxing  a  bug  when  it  is  encountered.  This  means  that  after  i 
faults  have  been  found,  we  expect  i  x  p  faults  to  have  been 
ctxrected,  instead  of  i.  Thus  the  failure  rate  of  T.  is 

r.j(tlNAp)  =  A(N-p(i-l)) 

i 

When  p  =1  we  get  the  Jelinski  &  Moranda  model. 

A  model  by  Schick  &  Wolverton  (1978). 

This  is  yet  another  Type  I  model,  and  this  time  the  failure 
rate  is  assumed  proportional  to  the  number  of  bugs  remaining 
in  the  system  and  the  time  elapsed  since  the  last  failure.  Thus 

r^(tlNA)  =  A(N-i+l)t,  t>0 

This  model  differs  from  models  M  in  that  the  failure  rate 
does  not  decrease  monotonically.  Immediately  after  the  i-th 
failure,  the  failure  rate  drops  to  0,  and  then  increases  linearly 
with  slqte  (N-i)  until  the  (i+l)th  failure;  see  Figure  2(iv). 

Bayesian  Differential  Debugging  Model 
(Littlewood  (1980)). 

This  model  can  be  considered  as  an  elaboration  of  model 
2  proposed  by  Littlewood  &  Verall.  Recall  that  in  model  2  it 
was  assumed  that  A.,  the  failure  rate  of  the  i-th  time  between 
failures,  was  declar^  to  have  a  gamma  distribution.  In  this 
new  model  Littlewood  supposed  that  there  were  N  bugs  in  the 
system  (a  return  to  the  bug  counting  phenomenon),  and  then 
proposed  that  A.  be  specified  as  a  function  of  the  remaining 
bugs.  In  particular,  he  slated  Aj  =  ^,  +  +  ...  +  where 

were  independent  and  identically  distributed  gamma  ran¬ 
dom  variables  with  shape  a  and  scale  3.  This  implied  that  A. 
would  have  a  gamma  distribution  with  shape  a(N-i)  and  scale 
3.  In  other  respects  its  assumptions  are  identical  to  the  original 
Littlewood/Veiall  model. 

Bayes  Empirical  Bayes  or  Hierarchical 
Model  (Mazzuchi  &  Soyer  (1988)). 


Compared  to  the  Jelinski  k  Moranda  model,  where  the 
(kop  in  failure  rale  after  each  figure  was  always  A,  the  drop  in 
failure  rale  here  idler  the  i-ih  failure  is  D  k''‘(l-k)  see  Figure 
2(iii).  The  assumption  of  a  perfect  fix,  with  no  introduction  of 
new  bugs  during  the  fix,  is  retained. 


In  1988  Mazzuchi  &  Soyer  proposed  a  Bayes  Empirical 
Bayes  or  Hierarchical  extension  to  the  Littlewood  &  Verall 
model.  As  with  the  original  model,  they  assumed  T  to  be 
exponentially  distributed  with  scale  A..  Then  they  proposed 
two  ideas  for  describing  A.,  here  called  model  A  and  model  B. 
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Model  A 

Still  assume  that  A.  is  described  by  a  gamma  distribution, 
but  with  parameters  a  and  p.  Now  assume  that  a  and  3  are 
independent  and  that  they  themselves  are  described  by  prob¬ 
ability  distributions;  a  by  a  unirorm  and  3  by  another  gamma. 
In  other  words: 


Ji  (a  I V)  =  ^,  0  ^  a  ^  V 


n(3la.b)  =  j^3->e-^,  3S0.a>0.b>0. 

Model  B 

Assume  that  A.  is  described  exactly  as  in  Littlewood  and 
Verall  i.e. 

n^(Xla,v(i))  =  ^^^X“-‘e-’^«\  X.>0 

and  that  y(i>s  3o  +  3i>>  except  now  place  probability  distribu¬ 
tions  on  01, 3o  and  3i  as  follows: 

n(a  I  (0)  =  — ,  0  ^  a  ^  0) 

w 

3„2-3,,a>0,b>0 

«(3,ic/i)=j^3,''‘c^^ 

3,20,c>0.d>0. 

So  a  is  described  by  a  uniform  distribution,  3^  by  a  shifted 
gamma  and  3t  by  another  gamma,  and  there  is  dependence 
between  3o  and  3]-  By  assuming  3,  to  be  degenerate  at  0, 
model  A  is  obtain^  from  model  B.  The  authors  were  able  to 
find  an  tqiproximation  to  the  expectation  of  T^.  given  that 
T,st,,  TjSlj.  -,  and  so  use  their  model  to  predict  future 

reliability  of  the  software  in  light  of  the  previous  failure  times. 

Time-dependent  Error  Detection  Model 
(Goel  &  Okumoto  (1979)). 

This  is  the  first  lype  II  model  that  we  will  consider.  It 
assumes  that  M(t),  the  number  of  failures  of  the  software  in 
time  [0,0,  is  described  by  a  Poisson  process  with  intensity 
function  given  by 

X  (t)  *  ab  e"** 


where  a  is  the  total  expected  number  of  bugs  in  the  system  and 
b  is  the  fault  detection  rate;  see  Figure  2(v).  Thus  the  expected 
number  of  failures  to  time  t  is: 

M.(t)=f  abe”‘*ds  =  a(l -e”**). 

0 

The  function  p(t)  completely  ^^ecifles  a  particular  Pois¬ 
son  process,  and  the  distribution  of  M(t)  is  given  by  the  well 
known  fmnula 

P  (M  (t)  =  n)  =  e~  W  n  =  0,1,2,... 

n! 

Experience  has  shown  that  often  the  rate  of  faults  in 
software  increases  initially  before  eventually  decreasing,  and 
so  in  Goel  (1983)  the  model  was  modified  to  account  for  this 
by  letting 

A.(t)  =  abct‘="’  e"**' 

where  a  is  still  the  total  number  of  bugs  and  b  and  c  describe 
the  quality  of  testing. 

Logarithmic  Poisson  Execution  Time 
Model  (Musa  and  Okumoto  (1984)). 

The  Logarithmic  Poisson  Execution  Time  Model  of  Musa 
and  Okumoto  is  one  of  the  more  popular  software  failure 
models  of  recent  years.  It  is  a  type  II  model,  but  the  model  is 
not  derived  by  directly  assuming  some  intensity  function  X(t), 
as  was  the  case  with  the  model  of  Goel  &  Okumuto.  Here  X(t) 
is  expressed  in  terms  of  p(t),  the  expected  number  of  failures 
in  time  [0,t)  via  the  relationship. 

A.(t)  =  Xoe-«'‘W. 

Put  simply,  this  relationship  encapsulates  the  belief  that 
the  intensity  (or  rate)  of  failures  at  time  t  decreases  exponen¬ 
tially  with  the  number  of  failures  experienced,  and  so  bugs 
fixed,  up  to  time  L  The  fixing  of  earlier  failures  will  reduce 
Ml)  more  than  the  fixing  of  later  ones.  Since  we  are  modeling 
the  number  of  failures  by  a  Poisson  process,  then  we  have 
another  relationship  between  X(t)  and  4(t),  namely 

H(t)  =  J  A.(s)ds. 

0 

Using  these  two  relationships  between  X,(t)  and  p(t),  there 
is  a  unique  solution  for  the  two  functions: 

Figure  2(vi)  shows  a  plot  of  X(t)  versus  t;  it  is  similar  to 
the  plot  of  figure  2(v)  except  that  the  tail  is  thicker. 

It  now  follows  from  the  above  that  by  using 
P  (M  (t)  =  n)  =  (p  (t))". ■  ^‘Vn!  we  can  say 
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P(M(t)  =  n)  = 


(ln(Xp9t+l))" 

e"(Xo0t+l)'^jcn!’ 


n  =  0,1.2... 


As  a  final  remark,  we  mention  that  in  their  ps^r  the 
authors  go  into  some  detail  on  estimation  of  and  6  by 
maximum  likelihood  methods;  however,  one  of  the  likelihoods 
aiqiears  to  be  incorrect. 

Random  Coefficient  Autoregressive 
Process  Model  (Singpurwalla  &  Soyer 
(1985)). 

This  is  a  IVpe  1-2  model,  that  is  one  that  does  not  consider 
the  failure  rate  of  times  between  failure.  Instead  it  assumes  that 
there  is  some  pattern  between  successive  failure  times  and  that 
this  pattern  can  be  described  by  a  functional  relationship  between 
them.  The  authcxs  declare  this  relationship  to  be  of  the  form 

=  ,  i  =  i,2,3,... 

where  To  is  the  time  to  the  first  failure  and  6i  is  some  unknown 
coefficient  If  all  the  6i’s  are  bigger  than  1  then  we  expect 
successive  lifelengths  to  increase,  and  if  all  the  6i’s  are  smaller 
than  1  we  expect  successive  lifelengths  to  decrease. 

Uncertainty  in  the  above  relationship  is  expressed  via  an 
error  term  5^,  so  that 


The  authors  then  make  the  following  assumptions,  which 
greatly  facilitate  the  analysis  of  this  model.  They  assume  the 
Tj’s  to  be  lognormally  distributed,  that  is  to  say  that  log  T.’s 
have  a  normal  distribution,  and  that  they  are  all  scaled  so  that 
T.  ^  1 .  The  5.’s  are  also  assumed  to  be  lognormal,  with  median 
1  and  variance  (the  conventional  notation  is  A(l,  a,^)). 
Then  by  taking  logs  on  the  relationship  above  they  obtain 

togT.  =  e.log7;_,+Iog& 

-“y- 


FHpuiw  3. 

Decision  process  for  single-stage  testing. 


Since  the  T.’s  and  the  6.’s  are  lognormal  so  the  log  T.’s 
and  the  e.’s  (=log  5.’s)  will  be  normally  distributed,  and  in 
particular  e.  has  mean  0  and  some  variance  (the  conven¬ 
tional  notation  is  A^(0,  o^).  The  log-lifelengths  therefore  form 
what  is  known  as  an  autoregressive  process  of  order  1  with 
random  coefficients  0..  There  is  an  extensive  literature  on  such 
processes  which  can  now  be  used  on  this  model. 

All  that  remains  to  do  is  to  specify  0.,  and  the  authws 
consider  sev^al  alternative  models.  For  example,  one  could 
make  0.  itself  an  autoregressive  process: 

0.  =  a0.  ,  +  ©.  ,  where  (O.isN  (O,  W.)  with  W.  known. 

When  a  is  known,  the  expressions  for  log  T.,  and  0. 
together  form  a  Kalman  filter  model,  on  which  there  is  also  an 
extensive  literature.  When  a  is  not  known  the  solution  is  via 
an  adaptive  Kalman  filter  algorithm  for  which  the  authors 
prqrose  an  approach.  As  an  alternative  to  the  above,  one  could 
place  a  two  stage  distribution  on  0. ,  and  the  authors  considered 
the  idea  of  0.  being  N(X,  with  X  also  a  normal  random 
variable  having  mean  m^  and  variance  s^^.  In  this  latter  case 
one  can  employ  standard  hierarchical  Bayesian  inference  tech¬ 
niques  to  i»edict  future  reliability  in  the  light  of  previous 
failure  data. 


4.  Model  Unification. 

By  adopting  a  Bayesian  approach,  it  turns  out  that  one  can 
unify  models  1, 2  and  8  -  the  models  by  Jelinski  &  Moranda, 
Littlewood  &  Verall  and  Goel  &  Okumoto  respectively  -  under 
a  general  framework.  Observe  that  this  also  provides  a  link 
between  the  two  types  of  models,  since  models  1  and  2  ate  of 
type  I  whilst  model  8  is  of  type  11. 

We  begin  by  recalling  the  first  model,  that  by  Jelinski  & 
Moranda.  Each  T.  is  assumed  to  have  a  constant  failure  rate 
A(N-i-«-l).  It  is  well  known  that  this  implies  each  T.  must 
therefore  be  exponentially  distributed  with  mean  (A(N-i-i- 1))  ’ . 
Now  assume  that  A  and/or  N  is  unknown;  in  true  Bayesian 
fashion  prior  disuibutions  are  placed  upon  them. 

To  obtain  model  8  by  Goel  &  Okomuto,  we  let  A  be 
degenerate  at  X  and  N  have  a  Poisson  distribution  with  mean 
0.  One  can  calculate  M(t)  using  the  T.’s  as  defined  by  Jelinski 
&  Moranda,  and  then  by  averaging  out  over  N  one  Ends  that 
M(t)  is  indeed  a  Poisson  process  with  mean: 

p(t)  =  0(l-e-^) 

which  is  the  form  of  p(t)  for  Goel  &  Okumoto ’s  model. 

One  can  also  obtain  model  2  by  assuming  N  to  be  degen¬ 
erate  and  to  have  a  gamma  distribution.  The  derivations 
which  lead  to  the  above  are  complex;  readers  are  referred  to 
Landberg  and  Singpurwalla  (198S)  for  the  details. 
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5.  An  Application:  Optimal 
Testing  of  Software. 

The  failure  models  that  have  been  reviewed  in  the  preced¬ 
ing  sections  can  be  used  for  more  than  inference  or  the  predic¬ 
tion  of  software  failure.  They  can  also  be  applied  in  the 
framework  of  decision  theory  to  solve  decision  {voblems.  An 
important  example  of  such  a  problem  is  the  optimal  time  to 
test  software  before  releasing  it  This  involves  the  balancing 
of  the  costs  of  testing  and  the  risk  of  software  obsolescence 
with  the  cost  of  in-service  failure,  should  a  bug  not  be  cor¬ 
rected  during  the  testing  period.  The  following  is  taken  from 
Singpurwalla  (1991),  in  which  a  strongly  Bayesian  approach 
is  taken. 

To  implement  a  decision  theoretic  procedure  requires  two 
key  ingredients.  The  first  is  a  probability  model,  and  here  we 
take  a  generalization  of  the  Jelinski  &  Moranda  model.  The 
second  is  a  consideration  of  the  costs  and  benefits,  or  utilities, 
associated  with  a  particular  decision  i.e.  the  costs  of  testing, 
the  benefits  and  costs  of  fixing  a  bug  etc.  Decision  theory  states 
that  the  optimal  decision  (in  this  case  time  of  test)  is  that  which 
maximizes  expected  utility. 

If  the  software  is  to  be  tested  for  some  time,  say  T  units, 
and  then  released  the  problem  is  to  find  a  T  that  maximizes 
expected  utility.  This  is  called  single  stage  testing.  There  is  a 
more  complex,  yet  realistic,  scenario  called  two  stage  testing. 
Here  the  software  is  tested  for  T  units  of  time,  and  then 
depending  on  how  many  failures  M(T)  were  observed  during 
that  test,  a  decision  is  made  on  whether  to  continue  testing  for 
a  further  T*  units.  The  problem  here  is  to  find  the  optimal  T 
and  T*,  with  T*  to  be  determined  before  M(T)  is  observed. 
Finally  there  is  a  third  testing  scenario,  namely  sequential 
testing.  Here  T*  is  determined  after  M(T)  is  observed;  this 
procedure  can  continue  for  several  stages,  with  T**  being 
determined  after  M(T*)  is  observed  and  so  on.  Here  we 
consider  the  case  of  single  stage  testing.  Figure  3  is  a  graph  of 
the  decision  process  associated  with  single  stage  testing. 

The  model  chosen  in  this  paper  is  an  extension  to  Jelinski 
&  Moranda ’s  model.  We  have 

(tlN,A)  =  A(N-i+ t>0 

In  the  previous  section  we  placed  prior  distributions  on 
one  of  N  or  A.  Now  we  place  priors  on  both  the  parameters, 
and  say  that  N  has  a  Pdisson  distribution  with  mean  9,  A  has 
a  gamma  distribution  with  scale  p  ^d  shape  a  and  that  N  and 
A  are  indqiendent. 

We  now  turn  to  the  choice  of  utility  function.  The  follow¬ 
ing  assumptions  are  made: 

i)  The  utility  of  a  program  that  encounters  j  bugs  during 
itt  operation  is  a,  +  a2  e~*i'. 

iO  The  cost  of  fixing  a  bug  is  some  constant  C, 
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Time  of  testing  versus  expected  utility  for  the  software  testing 
model. 


Expected  utility 


iii)  Let  f(T)  be  the  cost  of  testing  and  lost  c^poitunity  to 
time  t;  here  we  say  fCT)  =  dT*. 

Note  from  i)  that  the  utility  of  a  bug-free  program  is  a^+a^, 
and  the  utility  of  a  program  with  a  very  large  number  of  bugs 
is  near  a,,  so  that  typically  a,  is  a  large  negative  number 
(because  there  is  a  great  loss  associated  with  software  that  is 
constantly  failing  in  the  marketplace)  and  a,  >  0.  Combining 
these  assumptions  gives  us  the  utility  of  a  program  that  is 
tested  for  T  units  of  time,  during  which  M(T)  bugs  are  found 
and  ccxrected,  and  then  released  where)  bugs  are  encountered 
by  the  customer  as 

U  {T,M  (7) j)  =  e-*^  X  [a,  +  a.  e-T>  -  C,M  (7)  -  df] 

where  e'^^  is  some  devaluating  factor. 

Now  the  two  parts  of  the  decision  process  -  the  probability 
model  and  the  utility  function  -  are  brought  together.  We  wish 
to  find  the  time  f  that  maximizes  expected  utility.  In  other 
words  find  t  such  thatE(l/(T,  M(T),  j))  is  a  maximum,  where 
we  take  expectation,  using  our  failure  model,  with  respect  to 
M(T)  and  j.  This  maximization  is  quite  complex,  and  must  be 
done  numerically  via  computer.  The  details  are  found  in  the 
paper,  but  the  end  result  is  best  displayed  as  a  graph  of  time 
against  expected  utility  (figure  4);  in  this  case  one  can  sec  that 
the  time  one  should  test  the  software  for  is  about  3.S  units. 
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6.  Conclusion. 

This  paper  has  attempted  to  review  the  main  methods,  and 
some  of  the  more  well  known  models,  that  have  been  used  by 
the  statistics  community  in  the  area  of  software  reliability.  The 
Hrst  models  were  almost  always  based  on  looking  at  the  failure 
rate  of  the  software;  later  on  the  idea  of  modeling  number  of 
failures  by  a  Poisson  process  was  used  and  then  most  recently 
auto-regressive  processes  have  been  suggested  as  an  alterna¬ 
tive  to  the  failure  rate  method.  Application  of  the  failure 
models,  such  as  to  the  optimal  testing  decision  problem,  is 
another  important  aspect  to  the  field. 

Earner  it  was  pointed  out  that  there  is  almost  no  data  on 
the  reliability  of  commercial  software,  due  to  the  sensitive 
nature  of  that  information.  A  possible  method  of  overcoming 
this  problem  would  be  to  have  more  interaction  between  the 
statistics  and  computer  science  communities.  In  the  future, 
such  interaction  seems  essential  if  models  are  to  become  more 
realistic  and  useful,  and  it  is  perhaps  surprising  that  there  are 
so  few  links  between  the  two  groups  today. 

There  still  remains  much  to  be  researched  in  this  field.  In 
the  case  of  optimal  testing,  plans  for  two-stage  and  sequential 
testing  need  to  be  developed,  whilst  the  verification  of  current 
and  future  models  is  likely  to  remain  a  problem.  Nevertheless, 
because  of  the  increasing  presence  of  computers  in  all  aspects 
of  our  daily  lives,  the  tqiic  of  software  reliability  can  only 
become  more  impcHlant  in  the  future. 
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